Direct sampling of complex landscapes at low temperatures: 
the three-dimensional ± J Ising spin glass 
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A method is presented, which allows to sample directly low-temperature configurations of glassy 
systems, like spin glasses. The basic idea is to generate ground states and low lying excited configu- 
rations using a heuristic algorithm. Then, with the help of microcanonical Monte Carlo simulations, 
more configurations are found, clusters of configurations are determined and entropies evaluated. 
Finally equilibrium configuration are randomly sampled with proper Gibbs-Boltzmann weights. 

The method is applied to three-dimensional Ising spin glasses with ± J interactions and tempera- 
tures T < 0.5. The low-temperature behavior of this model is characterized by evaluating different 
overlap quantities, exhibiting a complex low-energy landscape for T > 0, while the T = behavior 
appears to be less complex. 

PACS numbers: PACS Numbers: 75.10.Nr, 75.50.Lk, 75.40.Mg, 05.50.+q 
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I. INTRODUCTION 

Despite large efforts made by the scientists in the last 
two decades, complex energy landscapes with many lpr 
cal minima and nested valleys, like that of spin glassesa, 
still offer many relevant questions to be answered. These 
questions usually regard the lowest energy levels of the 
landscape. The traditional numericaLapproach is to ap- 
ply a Monte Carlo (MC) simulation!]. Equilibration is 
tested by monitoring different average quantities as a 
function of the number of MC steps. Equilibration can 
be assumed, when the measured values of different runs, 
initially being Jjar apart, agree within error bars. An- 
other approaches is to calculate one quantity, like the link 
overlap, in two different ways, one time directly and one 
time depending on some other measured quantity like the 
energy, and wait till both results agree. 

Such a test is available only in special cases, e.g. for 
spin glasses with a Gaussian distribution of the bonds. 
Otherwise, one usually waits till the quantity of inter- 
est does not show any more a time dependence. Never- 
theless, at low temperatures and with increasing system 
size, equilibration becomes much harder and eventually, 
at very low temperatures, is impossible. 

In the very last years, a different approach has been 
proposed, namely the calculation of ground-state (GS) 
and low-energy configurations. Some characteristics of 
the low-energy landscape can be probed by the applica- 
tion of suitable perturbations which slightly modify the 
GSo. But the full information on the low-temperature be- 
havior can be obtained only by an equilibrium sampling 
of the system at a given temperature. Here we show, 
that by calculating GS and excited states, one can di- 
rectly sample rvery low temperatures. Several algorithms 
and heuristicsa are available to obtain ground states and 
excited states. Some are based again on Monte Carlo 
techniques like simulated annealing (SimA) and parallel 



tempering (PT). All these techniques have the drawback, 
that it is impossible to obtained an unbiased, i.e. equi- 
librium sample of configurations for T — » 0. For the MC 
methods, the reason is that for larger systems and very 
low temperatures, equilibration times are too long. We 
shall give below an example which shows for a ± J Ising 
spin glass, which exhibits an exponential ground state 
degeneracy, that just obtaining ground states is much 
easier than obtaining ground states with their proper 
statistics, i.e. each ground state with the same proba- 
bility. For other existing heuristics the statistics of the 
configurations is influenced in an uncontrollable way by 
the low-energy landscape. 

In this work, a post-processing method is presented, 
which removes the bias induced by the non-equilibrium 
low-temperature sampling and allows to obtain a prop- 
erly equilibrated state for systems having a high degener- 
acy. The basic idea of the technique is to calculate clus- 
ters of configurations, which arc connected in configura- 
tion space by zero-energy moves, e.g. zero-energy flips of 
spins in the Ising spin-glass case. Next, the sizes of these 
clusters are estimated and used to obtain an unbiased 
sample, where each cluster contributes with a factor to 
the size of the cluster and to the Gibbs-Boltzmann (G-B) 
weight. This method has already been successfully ap- 
plied to the ground-state sampling of three-dimensional 
Ising ±J spin glassesa. Here, the method is extended to 
the T > case and again applied to the d = 3 ±J SG 
model. Please note that this approach works better and 
better with decreasing temperature, hence is complemen- 
tary to the MC technique, which suffers from equilibra- 
tion problems at low temperatures. But similar to MC, 
one has to monitor some measured quantities as a func- 
tion of some parameters to establish equilibration, e.g. 
the number of clusters found in the analysis as a func- 
tion of the number of states included. Also similar to 
MC, obtaining equilibrium becomes harder with increas- 
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ing system size. In this sense, the method is also not ex- 
act. But in contrast to MC, ensuring equilibrium in this 
way is possible at very low temperatures for larger sys- 
tems (and becomes impossible for higher temperatures), 
while for MC it is the other way round. 

We apply the algorithm to three-dimensional Ising 
spin glasses. The EA model consists oi N — L 3 Ising 
spins Si = ±1 on a cubic lattice with the Hamiltonian 
H = — y^/j j\ JijSi-Sj. The sum runs over all pairs of 
nearest neighbors The Jij are quenched random 

variables taking values Jij = ±1 with equal probability 
and satisfy the constraint Y^,/ij) Jij — 0- We apply peri- 
odic boundary conditions in all directions. 

In this work we show that the overlap distribution 
P(q) at zero temperature is qualitatively different from 
P(q) at low but non-zero temperature. This means, even 
if there is an exponential number of GS configurations, 
zero-temperature quantities may be very different from 
those at any finite and small temperature. In particu- 
lar we will show here that for the three-dimensional EA 
model, which has a finite zero-temperature entropy, P(q) 
is very narrow at exactly T = 0, while it is broad at any 
finite temperature. We obtained the same result for the 
box-overlap Pbox(<?)- The picture resulting from our find- 
ings is that of a large number of GS which are however 
very close. Nevertheless, quite different states can be 
easily found once the first excited energy levels are con- 
sidered. This picture agrees with the very recent MC 
results by Palassini and Youngu. 

Before proceeding with our results and methods, we 
show, as a motivation, results from applying the SimA 
method to one sample realization of site L = 5 of our 
model. We have performed 10 4 independent runs of the 
SimA algorithm, starting with a temperature Tq = 2 and 
reducing the temperature according T n+ \ = bT n until 
T = 0.1 is reached. Per temperature 10 MC sweeps 
were performed. At the end of the simulation, one ran- 
domly chosen configuration exhibiting the lowest energy 
encountered during the run was stored. After having per- 
formed 10 runs, only the true ground states were kept. 
A GS configuration and its mirror image, obtained by re- 
versing all states, are treated as being equivalent. As it 
turns out, the system has 59 distinct GS configurations. 
In Fig. [j] histograms of the number of times each GS has 
been found are displayed for b = 0.5 and b = 0.99. One 
sees clearly that for b = 0.5 different GS configurations 
occur with different frequenciesa, i.e. not all appear with 
the same frequency as requested by the G-B distribution. 
When cooling much slower, i.e. with b = 0.99, all GS are 
almost equiprobable. This means that just finding GS 
configurations is much easier than finding each GS con- 
figuration with the correct probability. 

For system sizes just slightly larger than L = 5, the 
number of GS and excited states is already huge (e.g. 
~ 10 16 for L = 8). For this system sizes it is im- 
possible to obtain a histogram similar to the one pre- 
sented above. Consequently, it is impossible to determine 
whether all GS are sampled with the correct statistics. 
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FIG. 1: Histogram of the number of times each GS is found 
with a SimA simulation of 10 4 independent runs for one L — 5 
realization of a ± J Ising spin glass. The temperature was de- 
creased according to T n +i = bT n , with To = 2 until T = 0.1 
is reached. At each temperature 10 MC sweeps were per- 
formed. For the upper panel b = 0.5, while b = 0.99 for the 
lower panel. 



This is even more true for excited states. Please note that 
this is the same, for more elaborate algorithms like par- 
allel tempering!! Since, as already pointed out, at very 
low temperatures and for system sizes like L = 10 it is 
impossible to equilibrate the system, other methods have 
to be applied. In this paper, we present a post-processing 
tool, which allows to correct the bias imposed by any al- 
gorithm and leads to an equilibrated sample. For sizes 
up to L = 10 and low temperatures up to T < 0.5 the 
additional effort is moderate, because only the few lowest 
levels of excited states have to be considered. For larger 
temperatures, the post-processing methods becomes in- 
tractable, but then conventional MC methods can be eas- 
ily applied. 

The rest of the paper is organized as follows. First, 
we explain the algorithms we have applied. In the next 
section, we present the result for the three-dimensional 
±J spin glass. Finally, a summary and a discussion are 
given. 
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II. ALGORITHMS 

The technique to obtained an equilibrated low-tempe- 
rature sampling consists of four steps: 

1 . Generate configurations for GS and the lowest lev- 
els of excitations. 

2. On each energy level: group configurations into 
clusters. 

3. Calculate sizes of clusters. 

4. Generate a sample of states for given temperature 
T, where each cluster contributes with a weight 
proportional to its size and to the G-B factor 
exp(—E/T), where E is the energy of the config- 
urations in that cluster. 

Now all four steps are explained. 

The basic method used here to generate the configu- 
rations, is the cluster-exact approximation (CEA) -tech- 
niques, which is a discrete optimization methocH de- 
signed especially fa r . ap jin glasses. In combination with 
a genetic-,algorithmEJE3 this method is able to calculate 
true GSfU up to L = 14, as well as excited configurations 
as a byproduct. Since the CEA technique is well estab- 
lished and described in several sources, the details are 
skipped here. For each system and each energy level, we 
have generated 1000 configurations with the pure genetic 
CEA algorithms. We will show below that this number 
of configurations is sufficient up to L = 10 and T = 0.5. 

By applying pure genetic CEA, ope does not obtain 
the true thermodynamic distributiorJIfl, i.e. not all con- 
figurations with the same energy contribute to physical 
quantities proportional to the G-B weight. This means 
the genetic CEA algorithm is biased. For small system 
sizes up to L — 4 it is possible to avoid the problem by 
generating all low-energy configurations; averages can be 
performed simply by considering each configuration once, 
weighted with the G-B factor. Since the degeneracy in- 
creases exponentially with the number N of spins and 
grows also strongly with the energy level, a complete enu- 
meration is not possible for larger system sizes or higher 
energies. Instead, one has to choose a subset of all con- 
figurations, where each configurations contributes with 
a probability proportional to the G-B weight. The pro- 
cedure described here, consisting of steps 2-4 mentioned 
above, is applied to ensure that all configurations appear 
with the correct probability in this selection. Please note 
that the following methods works for any set of states, 
independently of the method which has been applied to 
generate the states. I.e. also the results of many inde- 
pendent runs of a low-temperature MC simulation can 
be treated, in case an equilibration was not possible, e.g. 
for very low temperatures and larger system sizes. 

In step 2 of our method, we group the configura- 
tions into, clusters by performing the ballistic-search al- 
gorithmic: All configurations which are accessible via 
flipping of spins having zero local field (called free spins 



in the following), i.e. without changing the energy E, are 
considered to be in the same cluster. Please note that 
the Hamiltonian is symmetrical with respect to flipping 
all spins simultaneously Hence, for the rest of the paper 
and for all analysis steps, a configuration and its mirror 
image are regarded as being identical. The final result 
is a list of different clusters whose sizes are estimated as 
explained below. This list does not change if more than 
one configuration was initially found in the same cluster, 
since these cases are recognized and correctly handled. 
For completeness and to convince the reader that the 
method indeed works, we present some details in the fol- 
lowing. 

The algorithm is applied independently for all config- 
urations having the same energy. The starting point is a 
set of ng configurations. For clarity, first the straight- 
forward method to obtain the cluster structure is ex- 
plained. This method will not be applied finally. Af- 
terward, the method actually used is exposed. 

The straight-forward construction starts with one arbi- 
trary configuration. It is the first member of the cluster. 
All configurations which differ only by the orientation of 
one free spin are called neighbors. All the neighbors of 
the starting configuration are added to the cluster. These 
neighbors are treated recursively in the same way: All 
their neighbors which are yet not included in the cluster 
are added, etc. After the construction of one cluster is 
completed the construction of the next one starts with a 
configuration, which has not been visited so far. 

The construction of the clusters needs only linear 
computer-time as function o£jis (0( n s))> similar to the 
Hoshen-Kopelman techniqualj, because each configura- 
tion is visited only once. Unfortunately the detection 
of all neighbors, which has to be performed at the be- 
ginning, is of O(rig) since all pairs of states have to be 
compared. Even worse, all existing configurations of a 
given energy must have been calculated before. As e.g. 
a 5 3 system may exhibit already more than 10 5 GS and 
much more excited states, this algorithm is not suitable. 

Instead we use the following technique, based on the 
ballistic-search algorithmic. The basic idea of ballistic 
search is to use a test, which tells whether two configura- 
tions are in the same cluster. The test works as follows: 
Given two independent replicas {erf} and {erf} let D 
be the set of spins, which are different in both states: 
D = {i\af ^ erf }. Now BS tries to build a path of suc- 
cessive flips of free spins, which leads from {erf} to {erf} 
while using only spins from D. In the simplest version it- 
eratively a free spin is selected randomly from D, flipped 
and removed from D. This test does not guarantee to 
find a path between two configurations which belong to 
the same cluster, since it may depend on the order the 
spins are selected whether a path is found or not. But, 
if a path is found, then it is sure that both configura- 
tions belong to the same cluster. On the other hand, 
if both configurations belong to the same cluster, then 
the method finds a path with a certain probability which 
depends on the size of D. It turns out that the proba- 



4 



bility decreases monotonically with \D\. For example for 
N = 8 3 the method finds a path in 90% of all cases if the 
two. states differ by 34 spins. More analysis can be found 

The algorithm for the identification of clusters utilizes 
a collective effect, to overcome the problem that some- 
times a path is not found, even if two configurations be- 
long to the same cluster. It works as follows: the basic 
idea is to let a configuration represent that part of a clus- 
ter which can be found using BS with a high probability 
by starting at this configuration. If a cluster is large it 
has to be represented by a collection of states, such that 
the whole cluster is "covered". For example a typical 
cluster of a 8 3 spin glass consisting of 10 16 ground states 
is usually represented by only some few ground states 
(e.g. two or three). A detailed analysis of how many 
representing configurations are needed-as a function of 
cluster and system size can be found inll3. The details of 
the algorithm are as follows: in memory a set of clusters 
consisting each of a set of representing configurations is 
stored. At the beginning the cluster set is empty. Iter- 
atively all available configurations {<7j} are treated: For 
all representing configurations the BS algorithm tries to 
find a path to the current configuration or to its inverse. 
If no path is found, a new cluster is created, which is 
represented by the actual configuration treated. If {at} 
is found to be in exactly one cluster nothing special hap- 
pens. If {<Ji} is found to be in more than one cluster, 
it is called a bridge configuration and all these clusters 
are merged into one single cluster, which is now repre- 
sented by the union of the states which have represented 
all clusters affected by the merge. After all configurations 
have been treated the whole process is run again with the 
obtained set of clusters. This allows to find bridge config- 
urations which have not identified in the first iteration, 
because accidentally only one cluster had been created 
during the jfixst iteration, at the time the configuration 
was treatedlij. 

The BS identification algorithm has some advantages 
in comparison with the straight-forward method: since 
each ground-state configuration represents many ground 
states, the method does not need to compare all pairs 
of states. Each state is compared only to a few num- 
ber of representing configurations. Thus, the computer 
time needed for thC|-ealculation grows only a little bit 
faster than 0{n§ncf^, where nc is the number of clus- 
ters, which is much smaller than us- Consequently, large 
sets of configurations, which appear already for small sys- 
tem sizes like N — 5 3 , can be treated. Furthermore, the 
cluster structure of even larger systems can be analyzed, 
since it is sufficient to calculate a small number of con- 
figurations per cluster. The main point is that one has 
to be sure that all clusters are identified correctly. This 
is not guaranteed immediately, since for two configura- 
tions belonging to the same cluster there is just a certain 
probability that a path of free flipping spins connect- 
ing them is found. But this poses no problem, because 
once at least one state of a cluster has been found, many 



more states can be obtained easily by just performing a 
E =const Monte-Carlo simulation starting with the ini- 
tial state. Hence, one can increasing the number of states 
available quickly. The probability that all clusters have 
been identified correctly approaches very quickly unity 
with increasingjipmber of available states. Detailed tests 
can be found iril3. For all results presented here, we have 
checked that the clusters do not change when doubling 
the number of states. 
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FIG. 2: Cluster-size distributions of GS clusters for small sizes 
L — 3 to L — 8. The straight line represents the function 
2V- 1A . 



Furthermore, one has in principle to ensure that really 
all clusters are found, which is simply done by calculat- 
ing enough configurations, hut this is still only a tiny 
fraction of all configurationseJ. This time, the configura- 
tions must be obtained independently, one cannot use the 
E =const MC simulation as above. It is possible to ob- 
tain at least one configuration from each cluster roughly 
up to size L = 8 at GS level, resp. L = 6 for first excited 
states. For sizes like N — 10 3 , the largest size we have 
treated in this paper, the number of clusters is too large 
at any energy level. But this is not a problem in principle 
because the low-temperature behavior of these systems is 
dominated by large clusters. As an example, in Fig. |^ 
the probability densities of cluster sizes for GS clusters 
are shown. The distributions are for small system sizes 
up to L — 8, were we-can be fairly sureEJ that all clus- 
ters have been foundEl The distributions follow roughly 
an algebraic decrease with a p(V) ~ V~ a behavior with 
a ~ 1.1. This dependence gets straighter with increasing 
system size. We are interested in the contribution of a 
cluster of order (or scale) of size V to the behavior. First, 
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the statistical weight of a cluster is proportional to the 
number of states in the cluster, i.e. to the volume V. Sec- 
ond, each scale of cluster sizes contributes proportional 
to the scale itself, because we are integrating over all clus- 
ters of a given scale, i.e. this weight is also proportional 
to V. (Or in other words, to translate the probability 
densities into probabilities on a logarithmic scale, one 
has to multiply with V.) In total, clusters of sizes with 
scale V contribute with weight V 2 p(V) = V 2 ~ a . Since 
a « 1.1 < 2, the largest scale clusters dominate the be- 
havior. On the other hand, since p(V) rapidly decreases, 
the number of these dominating clusters is rather small, 
i.e. it is rather simple to obtain an equilibrated sample of 
configurations. For the first excited level we have found 
a = 1.3 < 2, while at higher excited levels the number of 
clusters is too large to really find all of them. This results 
indicates that at higher levels the distribution becomes 
broader, which limits the application of the method to 
the lowest level of excitations. This effect is studied be- 
low with more detail. We have restricted our analysis to 
the first 4 levels of excited states. 

Please note that the CEA method generates configura- 
tions from larger clusters with larger probabilityca, hence 
the large and important clusters are encountered on av- 
erage first in the calculations. For the system sizes we 
have treated here, except L = 10 and T — 0.5, about 
90% of all contributing states are typically from the top 
5 largest clusters and further 5% from the next 5 largest 
clusters. Then with the 1000 configurations we generated 
per energy level, we encounter typically up to 100 clus- 
ters, and we can be pretty sure that all thermodynamic 
relevant contributions are considered within the level of 
accuracy given by our statistical fluctuations. Only the 
results for L = 10 and T = 0.5, where higher level exci- 
tations contribute significantly, may not be equilibrated. 
This is demonstrated at the end of this section, after we 
have presented the remaining parts of our algorithm. 

The third step in the algorithm is the estimation of the 
cluster sizes. This works as follows. Let C be a cluster 
we want to measure in size and let's consider a random 
'reference configuration' {r^} belonging to this cluster. 
We define a test Hamiltonian H[s] = — J^i r i s i f° r { s i} 6 
C, being E{(3) and S(J3) the average extensive energy and 
entropy at inverse temperature (3. Then the size of C is 
given by exp[5(0)]. Since the GS of this Hamiltonian 
is unique (it is the reference configuration), i.e. S(oo) — 
0, we obtain from the microcanonical definition of the 
temperature T = dE/dS 

f E(0) 

5(0) = S{0) - S(oo) = AS= /3dE = 

JE(oo) 

/>oo />oo 

= / [E - E{oo)} d/3 = / (E + N)df3 ,(1) 
Jo Jo 

where the previous last equality comes from an integra- 
tion by parts and the last equality from the substitu- 
tion E(oo) = —N. In order to calculate this integral, 
we actually perform a fast MC simulation restricted to 



configurations {s^ £ C while varying w = exp(— 2/3) 
in [0, 1] and measuring the average energy E as a func- 
tion of w. The final formula is the integral of a smooth 
function AS = / Q ^fdw. The number of MC sweeps 
applied per integration step was chosen automatically by 
the program in a way that the resulting entropy did not 
change by more than 5% of the value when the number of 
MC sweeps was doubled. I.e. the program started always 
with 10 MC sweeps, calculated the entropy integral, then 
applied 20 MC sweeps and so on. For small clusters, the 
calculation usually stopped after 20 MC sweeps. For the 
largest clusters encountered here, the algorithm stopped 
after the integration using 640 MC sweeps. We have also 
checked, that for these cases the measured entropy did 
not depend monotonically on the number of MC sweeps, 
i.e. we are sure that we did not miss a systematic trend 
when stopping the calculation at one point. 

In principle, there could be high entropic barriers, 
which prevent the size calculation from converging to the 
correct value. Fortunately, the full algorithm is not sus- 
ceptible to that problem. The reason is that the BS clus- 
tering method uses single spin flips at constant energy 
as well to determine the cluster structure, as described 
above. This means, if two parts of a cluster are con- 
nected through a very tiny path (the entropic barrier), 
which is not detected by the MC integration, the cluster- 
ing method is also not able to recognize both subclusters 
as belonging to the same cluster. Hence, if both sub- 
clusters are large, the genetic CEA method will have cal- 
culated with high probability configurations from both 
subclusters. In the analysis, because they are not identi- 
fied as belonging to the same cluster, they will appear as 
two independent large cluster, i.e. the correct statistics is 
ensured at the end. If on the other hand, one subcluster 
is small, it has a negligible contribution to the overall 
behavior, like other small clusters. 

After estimating the cluster sizes, a certain number of 
configurations is selected from each cluster, this is the 
last step of the algorithm listed in the beginning of this 
section. This number of configurations is proportional to 
the size of the cluster and to the G-B factor exp(—E/T). 
It means that each cluster contributes with its proper 
weight. This is possible for small temperatures and small 
sizes, where only few low-energy levels contribute to the 
thermodynamical behavior. 

The selection of the configurations is done in a manner 
that many small clusters may contribute as a collection 
as wellQ. For example, assume that 100 configurations 
are selected from a cluster consisting of 10 10 configura- 
tions, then for a set of 500 clusters of size 10 7 each (with 
the same energy) a total number of 50 configurations is 
selected, i.e. 0.1 configurations per cluster on average. 
The correct handling of such situations is achieved by 
first sorting all clusters in ascending order. Then the gen- 
eration of configurations starts with the smallest cluster. 
For each cluster the number of configurations generated 
is proportional to its size, to exp(— E/T) and to a fac- 
tor /. If the number of configurations grows too large, 
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only a certain fraction f% of the configurations which have 
already been selected is kept, the factor is recalculated 
(/ <— f * fi) and the process continues with the next 
cluster. 

The configurations representing the clusters are gen- 
erated from the initial configurations, obtained from the 
heuristic algorithm, by microcanonical MC simulation, 
i.e. iteratively spins are randomly selected and flipped if 
they are free. Since within a cluster there are no energy 
barriers, for the system sizes up to L = 10, applying 
100 MC sweeps ensures that all configurations within a 
cluster are visited with the same frequency. 

To summarize, by applying the algorithm presented 
here, each cluster appears with a weight proportional to 
its size and to exp(— E/T) and each configuration within 
a cluster appears with the same probability. Therefore, 
on total, the correct thermodynamic distribution is ob- 
tained. 




FIG. 3: Result for X0.5 (see Eq. (g) for definition) as a func- 
tion of the number N con { of configurations included per energy 
level in the analysis. The error bars at the right represent the 
limiting values N con f — > oo obtained from fitting the data 
points iV con f > 40 to algebraic functions. For small temper- 
atures T, few configurations are sufficient while at T = 0.5 
more than 1000 configurations are necessary. 



Nconf- The configurations were taken in the order they 
appeared in the generation using the genetic CEA, i.e. 
for a small number of configurations, the large clusters 
are more likely to be represented than the smaller clusters 
since genetic CEA preferentially generates configurations 
from larger clusters. One can see that for low temper- 
atures, even few generated configurations are sufficient 
to yield the true behavior. Please note that the remain- 
ing fluctuations are due to the fluctuations between the 
different samples of configurations. The reason that few 
configurations are sufficient here is that at low tempera- 
tures the GSs dominate and the number of GS clusters is 
fairly small. With increasing temperature, excited states 
become more important. For excited states, much more 
clusters exists. Thus, more configurations must be in- 
cluded into the analysis. This is visible in Fig. ||, where 
at e.g. T = 0.5 £0.5 depends strongly on -/Vc 0n f. For 
-^conf = 1000 T — 0.5 seems to be the borderline case, 
while for T < 0.5 the result for X0.5 seems to be converged 
(within error bars). We have checked this explicitly by 
fitting algebraic functions to the data points N con { > 40, 
resulting in an agreement within error bars of the limit- 
ing value N con { — > 00 with the result we have obtained at 
•Nconf - * 00. Hence we can be again confident that using 
1000 configurations per energy level, the results obtained 
here up to L = 10 and T < 0.5 represent the true equi- 
librium behavior or, at least, is so close to the true result 
that it cannot be distinguished from it at the level of 
accuracy determined by the statistical fluctuations. For 
smaller sizes, the number of clusters is smaller on each 
energy level, which means that 1000 configurations per 
realization and energy level are sufficient for even higher 
temperatures. But we restrict our analysis to T < 0.5 
here. 




energy level 



We have tested whether our generated data represents 
the equilibrium behavior by calculating the small-overlap 
weight X0.5, as defined in the beginning of the next sec- 
tion in Eq. (Q). £0.5 is obtained for the largest system 
size L = 10 and for different temperatures T as a func- 
tion of the number of configurations N con f included in the 
analysis per energy level. The result is shown in Figj^. 
Please note that the full analysis, as explained in this sec- 
tion, has to be repeated independently for each number 



FIG. 4: Fraction of configurations sampled from each energy 
level at T = 0.5 for different system sizes. Energy level is 
the ground state. Lines are guides to the eyes only. 



Finally, in Fig. ^, the fraction of configurations sam- 
pled at T = 0.5 for the different energy levels is shown for 
different system sizes. For the smallest size L = 4 almost 
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only GS configurations contribute to the thermodynam- 
ics, while increasing system size higher energy configu- 
rations become more important. Please note that only 
for L — 10 configurations from excitation level 3 con- 
tribute. There the degeneracy is much larger than for 
the lower levels. This explains, why the result for L = 10 
and T = 0.5 is probably not equilibrated. The result of 
Fig. |i| shows that, when studying the low-temperature 
behavior of glassy systems, it is not sufficient to study 
just GS configurations since the G-B factor and the size 
of the cluster (i.e. the entropy) must be taken into ac- 
count. Nevertheless for low temperatures and not too 
large system sizes, the energy levels which actually con- 
tribute to the partition function are very few. 



III. RESULTS 

We have calculated ground states and excited config- 
urations up to level four, for system sizes L < 10. Up 
to 3000 realizations of the disorder were considered (900 
for the largest system size). From the set of configura- 
tions, samples of several hundred equilibrium configura- 
tions were generated for temperatures T € [0, 0.5]. 



at finite temperatures. This can be seen even better, by 
calculating the fraction 




FIG. 5: Distribution P(\q\) of overlaps at T = 0.5 for different 
system sizes. Lines are guides to the eyes only. The inset 
shows the average weight X0.5 of the distribution for \q\ < 0.5 
as a function of system size for T = 0.5, 0.4, 0.3, 0.1. The lines 
represent fits to functions of the form x(L) — x°° + a L x , with 
x°° = and A = -1.10(5) for T = 0.1, x°° = 0.051(13) for 
T = 0.3, x°° = 0.095(4) for T = 0.4 and x°° = 0.122(4) for 
T = 0.5. 



For each disorder realization and each temperature, the 
distribution Pj(q) of overlaps q = i ^ s i* s f was calcu- 
lated, where {sf }, {s^} are two different equilibrium con- 
figurations. In Fig. H the disorder-averaged distribution 
p (k\) = [Pj(\<1\)]j is shown for T = 0.5, where [. . ]j de- 
notes the average over the quenched disorder. The long 
tail to q = seems to saturate at a finite weight, indi- 
cating the existence of a complex low-energy landscape 



P(q) dq 



(2) 



of overlaps smaller than qo. The result for qo — 0.5 is 
presented in the inset of Fig. ^| For zero temperatures, 
where only GS configurations-are sampled, xo.5 converges 
to or to a very small valueEl The rate of convergence 
is described by the finite-size dependence xq^,(V) ~ L x . 
We find A = —1.10(5), which is compatible with the_pre- 
dicted bound A < — 1 given by the "TNT"-scenarioM In 
Ref. |^ a larger value A = —0.90(10) was found. This 
slight difference might be due to the different ensembles 
studied, since in Refj^ the constraint Y^Uj) <Aj = was 
not applied. 

Please note that for small temperatures we sample only 
GS configurations, due to small system sizes. For larger 
temperatures T > 0.3, the asymptotic value of £0.5 is 
clearly larger than zero. Please note that the last point 
L = 10, T = 0.5 may not be converged, as discussed 
above. But, as you can see in Fig. ^, the value of X0.5 is 
an increasing function of the number of states included in 
the calculation. Hence, the true result (we have obtained 
Xo.f 10 (0.5) = 0.137(6) by extrapolating N con f — ► 00 as 
opposed to 0.126(7) found for N COD f = 1000) is probably 
above our value, thus supporting even more the conclu- 
sion that £0.5 > 0. 

Our results are quantitatively comparable to the data 
found in Ref. which were obtained by a parallel- 
tempering MC simulation. Although the authors had 
no reliable criterion to check equilibration of the system 
(in contrast to the case with Gaussian distribution of the 
disorderu), by comparison with our results it is very likely 
that in Ref. |?] indeed thermal equilibrium was obtained. 

A non-trivial distribution of overlaps is not a sufficient 
criterion for a complex energy landscape. A qualita- 
tively similar overlap distribution with a nonzero weight 
for small values of q would be obtained also for a sys- 
tem, where various configurations differ by a domain wall 
through the system at different positions, e.g. a ferromag- 
net .with antiperiodic boundary conditions in one direc- 
tionB. 

To rule out this scenario, we have calculated also the 
distributions of box (or window) overlapsB3IHJ. This over- 
lap is defined as usual, but restricted to a finite "window" 
of volume I x I x I, with I < L fixed independently of the 
system size L. Please note that for the aforementioned 
ferromagnet, the distribution of box overlaps converges 
to a pair of delta functions at q — ±1 when L — > 00. 
The result for I = 3, T = 0.5 is exhibited in Fig. || At 
finite temperature, similar to the conventional overlap, 
the low-g tails seems to saturate, but more slowly, at a 
non-zero weight with increasing systems size. This can 
be seen from the inset of Fig. ^, where X0.5 is shown as a 
function of system size for T = 0.1, 0.3, 0.4 and T = 0.5. 
For T > 0.3, X0.5 clearly converges to a nonzero value. 
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\q\ 

FIG. 6: Distribution -Pbox(M) of box overlaps at T = 0.5 
for different system sizes. Lines are guides to the eyes only. 
The inset shows the average weight xo.s of the distribution for 
\q\ < 0.5 as a function of system size for T = 0.5, 0.4, 0.3, 0.1. 
The lines represent fits to functions of the form x(L) = x^ ox + 
a b L x », with :r^, x = and X b = -0.86(5) for T = 0.1, x^ ox = 
0.05(13) for T = 0.3, a;^ x = 0.10(1) for T = 0.4 and a;£, x = 
0.13(1) for T = 0.5. 



Pi{qi) of link overlaps qi = s f s< j s i ' s j ■ Tne re ~ 

suit for T — 0.5 and different system sizes can be ob- 
served in Fig. [7| The distribution becomes narrower, 
but a second small peak seems to emerge. In the in- 
set of Fig. ^ the finite-size dependence of the variance 
(T 2 = f Q (q — q) 2 Pi(q) dq is shown for different tempera- 
tures. In all cases, the width seems to converge toward 
zero. Please note, however, that we cannot exclude that 
the variance converges to a small but finite value. When 
we fit it to a function of the form cr 2 (L) = a 2 ^ + a a L x ° 
we obtain, for T = 0.5, cr^ = 0.0038(28) with X 2 per 
degree of freedom of 0.1, which is a very good fit. Nev- 
ertheless, a Pi(qi) consisting of two peaks at distance of 
0.1 with weights 0.1 and 0.9 respectively has a variance 
cr 2 = 0.0009. 

The behavior of Pi{qi) is quantitatively the same for 
three-dimensional spin glasses with a Gaussian distri- 
bution of the interactional, which were found with a 
parallel-tempering MC simulation. 



Thus, we can conclude that indeed at finite tempera- 
tures, three-dimensional spin glasses exhibit a complex 
low-energy landscape. 

Please note that the non-trivial behavior occurs for 
low temperatures, probably for all temperatures T > 0, 
which are sufficiently far away from the phase transition 
T c w 1.1. Hence, the effects which were-found within a 
Migdal-Kadanoff approximation schemeEj are unlikely to 
explain the kind of behavior we find. 




0.2 0.4 0.6 0.8 1 
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FIG. 7: Distribution Pi(qi) of link overlaps at T = 0.5 for 
different system sizes L. Lines are guides to the eyes only. 
The inset shows the variance a 2 as a function of system size 
for T = 0.1,0.5. The lines represent fits to functions of the 
form a 2 {L) = a t L x > (L > 4), with A; = 0.53(5) for T = 0.1 
and A; = 0.27(1) for T = 0.5. 



Finally, we have computed the average distribution 



IV. SUMMARY 

Summarizing, we have presented an algorithm which 
allows to investigate the low-temperature behavior of 
Ising systems with high degeneracy by direct sampling of 
GS and excited configurations. The basic idea is to gen- 
erate configurations with any suitable algorithm, group 
the configurations into clusters, measure the size of the 
clusters and then obtain a very good estimate of the G-B 
measure, to sample configurations with. Similar to MC, 
where one has to increase the number of MC sweeps until 
the system is equilibrated, one has to increase the num- 
ber of independent configurations until the true behavior 
is obtained. The main difference to MC techniques is 
that the method presented here works better with de- 
creasing temperature, while MC equilibrates faster with 
increasing temperatures. In this sense these methods are 
complementary. 

We have applied the algorithm to study the low- 
temperature behavior of three-dimensional ±J Ising 
spin glasses. We find that the statistical properties 
of the exponentially many ground state configurations 
are not representative of the low-temperature behavior. 
In particular we have shown for the three-dimensional 
Edwards- Anderson model that both the distributions of 
the overlap and of the box-overlap seem to be very nar- 
row functions at T = 0, where only few states contribute 
to the G-B measure, and broad for finite T. Hence the 
model does have a complex state space, which seems to 
become trivial at T = 0. For this reason one is forced to 
probe the energy landscape at T > 0. The distribution 
of the link-overlap seems to develop a second peak, but 
the extrapolation of the asymptotic shape is beyond our 
present computational capabilities. 
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